Discrete breathers in nonlinear network models of proteins 
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We introduce a topology-based nonlinear network model of protein dynamics with the aim of 
investigating the interplay of spatial disorder and nonlinearity. We show that spontaneous local- 
ization of energy occurs generically and is a site-dependent process. Localized modes of nonlinear 
origin form spontaneously in the stiffest parts of the structure and display site-dependent activa- 
tion energies. Our results provide a straightforward way for understanding the recently discovered 
link between protein local stiffness and enzymatic activity. They strongly suggest that nonlinear 
phenomena may play an important role in enzyme function, allowing for energy storage during the 
catalytic process. 
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The predictions of clastic network models (ENMs) of 
proteins [TJ [2 |3J H] have proven useful in quantita- 
tively describing amino-acid fluctuations at room tem- 
perature [TJ, often in good agreement with isotropic [2J, 
as well as anisotropic measurements [5J [B]- Moreover, 
it has been shown that a few low-frequency normal 
modes can provide fair insight on the large amplitude 
motions of proteins upon ligand binding [71 [51 [S] , as pre- 
viously noticed when more detailed models were consid- 
ered [lOj [TTJ [12] , also by virtue of the robust character of 
the collective functional motions [13 . 

However, low-frequency modes of proteins are known 
to be highly anharmonic [141 115] , a property which has to 
be taken into account in order to understand energy stor- 
age and transfer within their structure as a consequence 
of ligand binding, chemical reaction, etc [TBJIT7]. Indeed, 
there is growing experimental evidence that long-lived 
modes of nonlinear origin may exist in proteins [T51 119j . 
Likewise, many theoretical studies have appeared sug- 
gesting that localized vibrations may play an active role 
in, e.g., enzyme catalysis [2Uj. These include topolog- 
ical excitations such as solitons [3JJ as well as discrete 
breathers (DBs) [22j[23]. 

The latter are nonlinear modes that emerge in many 
contexts as a result of both nonlinearity and discrete- 
ness Although their existence and stability proper- 
ties are well understood in systems with translational in- 
variancc, much less is known of the subtle effects arising 
from the interplay of spatial disorder and anharmonic- 
ity [551 HE1 127] ■ For this purpose, in the present work 
we introduce the nonlinear network model (NNM). Our 
aim is to extend the simple scheme of ENMs, known to 
capture the topology-based features of protein dynam- 
ics [TJ [51 [3J, by adding anharmonic terms. Within the 
NNM framework, we show that spontaneous localization 
of energy can occur in protein-like systems and that its 
properties may be intuitively rationalized in the context 
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Figure 1: Energy as a function of time, when citrate synthase 
is cooled down as a consequence of surface friction. Dashed 
line: total energy. Solid line: energy of Threonine 208, the 
amino-acid the most involved in the DB. Dotted line: energy 
of Alanine 209, also involved in the DB. k B T eq = 20 kcal/mol. 

of specific biological functions. In our model, the poten- 
tial energy of a protein, E pi has the following form: 
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where dij is the distance between atoms i and j, d®j 
their distance in the structure under examination (as e.g. 
solved through X-ray crystallography) and R c is a cutoff 
that specifies the interacting pairs. As done in numerous 
studies, only C a atoms are taken into account [3] and 
k% is assumed to be the same for all interacting atom 
pairs [TJ. As in previous ENM studies [HI HE]> we take 
R c =10 A, and fix k,2 so that the low-frequency part of 
the linear spectrum match actual protein frequencies, as 
calculated through realistic force fields pjjj HU [T5 . This 
gives &2 = 5 kcal/mol/ A 2 , with the mass of each C a 



2 




40 60 80 100 

Frequency (cm-1) 

Figure 2: Locality of citrate synthase harmonic modes, as a 
function of their frequencies, together with the locality and 
frequency of a discrete breather (DB). 

fixed to 110 a.m.u., that is, the average mass of amino- 
acid residues. Note that standard ENM corresponds to 
fej = 0, while in the present work k\ = 5 kcal/mol/A 4 . 

Proteins live and perform their functions immersed in 
water and exchange energy with the solvent through their 
sizable surface portion. In a previous paper we showed 
that complex energy relaxation patterns are observed as 
a result of the inhomogeneity of the coupling to the sol- 
vent of bulk and surface atoms In the presence 
of nonlinearity, boundary relaxation is known to drive 
a wide array of systems towards regions of phase space 
corresponding to localized modes that emerge sponta- 
neously [301 I3U 1321 133]- Thus, in order to study typical 
excitations of nonlinear origin in protein structures, it 
appears natural to perform a boundary cooling experi- 
ment. Our protocol is the following. After 50 psec of mi- 
crocanonical molecular dynamics (MD) simulation per- 
formed at a temperature T eq , the protein is cooled down 
by adding a linear dissipation term to the force acting 
on surface atoms, that is, those belonging to amino-acids 
with more than 25 A 2 of solvent accessible surface area. 
This represents nearly 40% of the amino-acid residues, 
for all proteins considered in the present study. The 
viscous friction coefficient 7 is set to 2 psec -1 , a typi- 
cal value for protein atoms in a water environment |16j . 
Hereafter, the equilibration energies considered are in the 
range k B T eq = 2 — 20 kcal/mol, that is, of the order of, 
e.g., the energy release of ATP hydrolysis. With such ini- 
tial conditions, energy in the system remains high for a 
period of time long enough so that localization can occur. 

In Fig. [T] we show the energy of dimeric citrate syn- 
thase (PDB code 1IXE) as a function of time, as well as 
the energy of two amino-acids of monomer A, Thr 208 
and Ala 209. After t = 20 psec and a few large fluc- 
tuations, a DB centered at Thr 208 forms. At t — 200 
psec, more than 80% of the total energy is located there. 
Note the slow decay of the total energy after t —100 psec 
and the periodic energy exchanges of Thr 208 with Ala 
209, another among the few amino-acids involved in the 
DB. Note also that at t = 20 psec the energy of Thr 
208 is higher than at t =0, that is, when the friction 
was turned on, a clear-cut demonstration of the known 
tendency of DBs to harvest energy from lower-energy ex- 
citations |24j . In order to check that the phenomenon 
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Figure 3: Stiffness of dimeric citrate synthase as a function 
of residue number (dashed line). The number of DBs found 
at a given site out of 500 instances is also reported (black 
diamonds, right y-axis). 



shown in Fig.JTJis indeed the spontaneous localization of 
a DB, we switched off the friction at t — 200 psec and 
performed 100 more psec of microcanonical MD simula- 
tion. Then, we projected the latter trajectory on the first 
eigenvector of the corresponding velocity-covariance ma- 
trix, which gives the pattern of correlated atomic veloci- 
ties involved in the DB. The Fourier transform of such a 
projection yields an accurate value for the DB frequency, 
while the spectral line-width provides information on the 
DB stability over the 100 psec analysis time-span. 

In Fig. [2] we report the harmonic spectrum of cit- 
rate synthase as well as the DB frequency as func- 
tions of a locality measure. The latter is defined as 

Lk =Ei,«fc*a] 4 /[Ei 1 a[e£J a ] a ' wh ere £f Q is the a (£,2,, z) 
coordinate of the i-th. atom in the fc-th displacement pat- 
tern (normalized eigenvector, DB). As expected, the DB 
frequency (130 cm -1 ) lies above the highest frequency of 
the harmonic spectrum (101 cm -1 ). Moreover, the cor- 
responding spatial pattern is much more localized than 
any of the harmonic modes (note the logarithmic scale). 

Starting from random initial conditions, we obtained 
500 stable DBs following the above-outlined protocol. Al- 
though in many cases several DBs emerged, we decided 
to retain only the runs where a single DB catched most 
of the system energy, and more energy than the aver- 
age amount per site at t = 0. In Fig. [3] we report the 
number of DBs found at each site. The largest fraction 
(20 %) of these highly energetic DBs formed at Thr 208 
in monomer A, but we also observed DBs at 27 other 
sites, noteworthy at Thr 192 of monomer B (18%). Note 
also that, although the studied protein is a dimer, that 
is, with a an approximate but clear two-fold structural 
symmetry, the probability to observe a DB at a given site 
varies from one monomer to the other, indicating that the 
localization dynamics is rather sensitive to small changes 
in the local environment. As shown in Fig. [3] this proba- 
bility is higher in the stiffest parts of the protein scaffold, 
as measured through an indicator of local stiffness s,. For 
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amino- acid i, the latter is denned as: 
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where Mi = J2j 6(R C — d%j) is the number of neighbors 
of the i-th residue and 6(x) is the Heaviside step func- 
tion. The second sum is over the set S of the ten high- 
est frequency harmonic modes. The averaging over the 
Mi neighbors slightly smoothes mode contributions and 
helps underlining the fact that in each monomer of cit- 
rate synthase there is a stretch of nearly fourty consec- 
utive amino-acids (residues 185-225) with a remarkably 
stiff environment, deeply buried at the interface between 
the two monomers. This is obviously where most DBs 
tend to emerge. Note, however, that the relationship be- 
tween high-frequency harmonic modes and spontaneous 
energy localization is not a straightforward one: for in- 
stance, DBs were observed only a couple of times at the 
site the most involved in the highest frequency normal 
mode, namely, Ser 213. As a matter of fact, as suggested 
by the large energy fluctuations observed at site Thr 208 
before the DB shown in Fig. [I] springs up, a competition 
between •potential DBs is likely to occur, with possible 
weak-to-strong energy transfers, before a given site is oc- 
cupied by a stable mode. 

In lattice systems sites are obviously equivalent. Here, 
as shown in Fig. |4j the energy-frequency relationship is 
site-dependent. Furthermore, the probability for a DB 
to localize at a given site depends in a non-obvious fash- 
ion upon the energy it needs to reach a given frequency 
at that location. While most DBs emerge at Thr 208, 
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Figure 4: DB frequencies in citrate synthase as a function 
of their energy (pluses). The cases of Threonine 208 (filled 
circles) and Alanine 209 (filled squares) are highlighted. Using 
our protocol, no DB with an energy lower than 37.8 kcal/mole 
was observed, out of a total of 500 cases. 
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Figure 5: Stiffness of the environment of amino-acid residues 
involved in enzymatic activity (black squares), compared to 
that of amino-acids of same chemical type (crosses) randomly 
chosen within the same set of enzyme structures. The broken 
line is only a guide for the eye. 

i.e. the site where the least energy is required for a given 
frequency, many DBs are also observed at Ala 209 in 
monomer A, one of the sites that demands more energy. 
In more than one dimension one expects DBs to appear 
only above a characteristic energy [2H EH]- Hence, our 
results hint at a strong site-dependence of such energy 
threshold, non-trivially related to local structural prop- 
erties. To shed light on this intriguing feature, a de- 
tailed characterization of the small-ampitude side of the 
DB energy-amplitude curves at different sites is currently 
under way. 

In the following step, we looked for DBs in other pro- 
teins, both dimeric and monomelic. For small proteins 
like HIV-1 protease (PDB code 1A30), a dimeric 2 x 99 
amino-acids enzyme, no DB could be obtained. This is 
likely to be due to the fact that in small proteins too 
many amino-acids are in direct interaction with a site 
where energy dissipation occurs. This means that small 
proteins may require more detailed models, like all-atom 
schemes, where cutoff values of the order of 5 A are 
customary [51 IT31 13"5] . 

In the case of aconitase (PDB code 1FGH), a 
monomeric 753 amino-acids enzyme, and alkaline phos- 
phatase (PDB code 1ALK), a dimeric 2 x 499 amino- 
acids enzyme, DBs prove nearly as easy to generate 
than in the case of citrate synthase. However, for pro- 
teins of similar sizes, the probability of similar events 
turns out to vary significantly from a protein to an- 
other. For instance, in the cases of phospholipase D 
(PDB code 1V0Y), a monomeric 504 amino-acids en- 
zyme, and isoamylase (PDB code 1BF2), a monomeric 
750 amino-acids enzyme, out of 100 cooling MD simu- 
lations, only 8 and 5 DBs were obtained, respectively, 
in contrast to citrate-synthase, where our success rate 
is over 50 %. This points to the intriguing conclusion 
that not only DBs in proteins are site-selective, but also 
appear to be non-trivially fold-selective. 
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In all the analyzed structures, spontaneous localiza- 
tion of energy occurs in the stiffest parts of the struc- 
ture. Thus, we turn now to examine the relationship 
between protein stiffness and function. Following the hy- 
potheses that enzymatic activity may require some kind 
of energy storage and that DBs may play a role in the 
process, we computed high-frequency normal modes for 
a set of 833 enzymes from the 2.1.11 version of the Cat- 
alytic Site Atlas Then, we determined the stiffness 
of each amino-acid known to be involved in enzymatic 
activity according to As a comparison, we also de- 
termined stiffnesses of amino-acids of the same chemical 
type, but picked at random among those not known to 
be involved in enzymatic activity. As shown in Fig. [5] 
catalytic amino-acids tend to be located in stiffer parts 
of enzyme structures, in agreement with our hypothe- 
ses. This is not an obvious result, since for the sake of 
catalytic activity amino-acids have to interact with en- 
zyme substrates, that is, to be accessible to them. Such 



a trend has already been noticed in other studies. Note- 
worthy, using the ease of displacing any given amino-acid 
residue with respect to the others as a stiffness measure, 
it was shown that roughly 80 % of the catalytic residues 
are located in stiff parts of enzyme structures |37j . In 
a more indirect way, it was also remarked that global 
hinge centers colocalize with catalytic sites in more than 
70 % of enzymes [3H] . So, stiff parts may play a role of 
pivot, allowing for accurate large-amplitude conforma- 
tional changes of enzymes upon substrate binding. 

What our results further suggest is that stiff parts of 
enzyme structures may also play another major role in 
enzyme function, namely by allowing for an active role 
of nonlinear localized modes in energy storage during the 
catalytic process. 

Y-.H.S. wishes to thank M. Peyrard and T. Dauxois 
for an invitation to talk at a training school held in Lcs 
Houches [25], where he was introduced to the fascinating 
world of discrete breathers. 
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